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Chapter 1 
INTRODUCTION 

1.1 Motivation 

A problem of considerable and increasing importance within the 
field of composite materials is the fracture of laminated composite 
structures. The anisotropy of the material greatly complicates even the 
simplist of problems. An example of a complicated problem is the edge 
replica of Fig. (1.1); it demonstrates the crack types present on the 
free edge of a [+302/-302] s tensile coupon. From Fig. (1.1) it is 
apparent that a crack can start out in a transverse mode and turn into a 
delamination within its growing length. Complexities such as this 
require in-depth fracture mechanics models which not only predict at 
what load a crack will extend but also the direction of crack exten- 
sion. This study was undertaken in an attempt to develop a model with 
the capability of describing the characteristics of crack growth in 
composites. 

1.2 Literature Review 

Smith [1] has discussed limitations of some of the current analyti- 
cal models for predicting crack growth characteristics in composite 
materials. Likewise, most of the previous finite element models either 
involve complex computational procedures or suffer from serious limita- 
tions. Some models distinguish between fiber and matrix [2, 3], others 
use hybrid or singular finite elements [4, 5], and still others assume a 
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direction of crack growth [6]. Finite element analyses which require a 
distinction between the fiber and matrix materials are impractical for 
most real life situations since the number of degrees of freedom 
required in such an analysis could easily exceed the capacity of most 
computers currently available. While the use of singular or hybrid 
finite elements may give a good representation of the stresses, strains 
and displacements near a crack tip, they generally require more computa- 
tion time and they are not well suited to explaining large scale crack 
growth characteri sties. Experimental results [7, 8] have shown that the 
direction of crack extension can change significantly during crack 
growth in a laminated composite material. These results limit analyses 
which assume a direction of crack growth to specialized cases or to 
small increments of crack growth. 

There are many different finite element models that do not require 
specialized elements. Two such models are the failed element approach 
of Adams [3], and the modified crack closure approach of Rybicki and 
Kannien [9]. The failed element approach assumes that when an element 
in an area of high stress exhausts its strain energy capacity, it 
fails. From this, it is assumed that a "crack" has formed and has the 
dimensions of the failed element. This approach has two implications, 
the most important of which is that a finite amount of material is 
removed from the system, which in an actual material is not the case. 

The other is that the crack is not likely to close up on itself in 
subsequent loading because of its exaggerated width. The modified crack 
closure technique is based on the crack closure integral and can be used 
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within the framework of a linear elastic analysis with a relatively 
coarse mesh. In the modified crack closure technique, the crack closure 
integral is evaluated directly from the nodal forces and displacements 
required to close a virtual crack of extension, Aa. The modified crack 
closure technique also has the advantage of obtaining mode-I, mode-II, 
and mode-III results in a single analysis. 

Few fracture theories predict the direction of crack extension [2], 
as well as the external load level which causes crack extension. Hashin 
[10] suggested that a failure criterion could be constructed which would 
include the plane on which failure would occur. Some of the many frac- 
ture/failure criteria which have been used to predict fracture/failure 
of composite materials include the Sih strain energy density criterion 
[2], the Tsai-Wu failure criterion [11], the Whitney-Nui smer point 
stress criterion [12], and the Hashin failure criterion [10]. Of the 
criteria listed, only the strain energy density criterion [2] and the 
point stress criterion [12] are readily capable of predicting the direc- 
tion of crack extension. However, without modification, their use is 
limited to special cases. 

1.3 Purpose of the Present Study 

The purpose of this study was to develop a finite element model 
capable of predicting crack growth characteristics in composite materi- 
als. It was desired to develop a model which not only could determine 
what applied load level would cause crack extension but one which could 
also determine the direction of crack extension. 
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1.4 Basic Assumptions 

Unless otherwise stated, the model developed was based on the 
following assumptions: 

(i ) Linear elastic, homogeneous isotropic or 
homogenous orthotropic fibrous composites 

(ii) small displacement theory 

(iii) a crack extends from one end only - one crack at 
a time 

(iv) no variation of geometry in one of the coordinate 
directions (i.e., plane stress, plane strain and 
generalized plane strain problems). 

1.5 Description of the Finite Element Model 

The finite element model developed uses a two dimensional mesh of 
four node, linear, isoparametric elements. The model has the capabili- 
ties of obtaining either plane strain, plane stress or generalized plane 
strain solutions. (Refer to Appendix A for a description of the finite 
element method as it applies to this study.) The material models avail- 
able include isotropic, orthotropic and laminated orthotropic, (off- 
axis), materials. (See Appendix B for an explanation of the constitu- 
tive relations for the respective material models.) 

The approach taken to the solution of the crack problem was to 
separate the analysis into two main parts. In the first part of the 
solution the direction of crack growth was determined and in the second 
part the load level which would cause crack extension was determined. 
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The crack growth direction was determined through the use of sev- 
eral fracture/fai lure theories. The theories considered include a 
modified version of the Griffith criterion [13], the Si h strain energy 
density criterion [2], the Tsai-Wu failure criterion [11], and modified 
versions of the Whitney and Nuismer point stress theory [12], and Hashin 
failure criteria [10]. The failure criteria are described in Chapter 2 
and the results are compared with theory and experiment in Chapter 3. 

The determination of the load level which would cause crack exten- 
sion was made through the use of the modified crack closure method 
[9]. The modified crack closure method, as it applied to this study, is 
presented in Chapter 2. 

The main reasons behind the choice of this solution approach was 
that it could be used with a linear elastic analysis, that it could be 
used with a relatively coarse mesh and that it required a minimum of 
computer time. 

1.6 Problems Considered 

The problems considered in this study were: 

(i) A mode-I crack in an infinite isotropic plate with remote 
loading of a. Fig. (1.2a), was considered as a test of 
the crack closure technique. 

(ii) Mixed mode cracks in infinite plates of isotropic materials. 

Fig. (1.2b), were analyzed as a test of the crack growth 

direction, -6 , for various angles of crack inclination, B. 
o 
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(iii) Off-axis, unidirectional, orthotropic tensile specimens. 

Fig. (1.3a), were analyzed for the direction of crack 
extension in fibrous composites and the results were compared 
with available experimental results. 

(iv) Transverse cracks were introduced on the free edge of a 

laminated composite tension specimen. Fig. (1.3b), and the 
predicted crack paths were compared against available 
experimental results. 
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a. Crack In A Unidirectional 
Tensile Specimen 



b. Transverse Crack In A Laminated 
Tensile Specimen 


Figure 1.3 Composite Crack Problems 
Considered 


Chapter 2 

THEORETICAL BACKGROUND 


2.1 Criterion for Crack Extension 

2.1.1 Energy Release Rate Concept 

In the Griffith theory of classical fracture mechanics [13], it is 
assumed that strain energy is released when a crack surface is created 
in a stressed body. The rate of energy release when a crack extends 
stably in a body is known as the critical energy release rate, G c> The 
critical energy release rate, G c , can be determined experimentally by a 
procedure which allows for stable, slow crack extension, (see Ref. [17] 
for isotropic materials and [6] for composites). 

For structures, such as a composite laminate, an existing crack may 
or may not grow under a given state of stress. To determine whether or 
not an existing crack will extend, it is necessary to calculate the 
available energy release rate, G(a), associated with a crack of length 
a. If the available energy release rate, G(a), is equal to the critical 
energy release rate, G c the crack will grow in a stable fashion. If the 
available energy release rate, G(a), is greater than the critical energy 
release rate, G c , the crack grows unstably and if the available energy 
release rate, G(a), is less than the critical energy release rate, G c , 
the crack does not extend. Similarly, the external load which first 
causes the available energy release rate, G(a), to reach the critical 
energy release rate value, G c , is the critical load and loads greater 
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than the critical load level result in unstable growth while loads less 
than the critical load level do not extend the crack. 

Mathematically, the energy release rate, G(a), for a given crack of 
initial length, a, is defined as the difference in the total strain 
energy of the structure, AU, before and after a small crack extension, 
Aa, is introduced, that is. 


G(a) 


lim AU 
Aa+0 Aa 


( 2 . 1 ) 


2.1.2 Modified Crack Closure Approach 

Irwin [14] contended that if a crack extends by a small 
amount, Aa, the energy released in the process is equal to the work 
required to close the crack back to its original length. This statement 
in equation form is 


G(a) = Aa™0 2AT / da (Z ’ 2) 

■f 

where a is the surface traction vector and Au the displacement vector 
required to close the crack back to its original length. 

The modified crack closure technique of Rybicki and Kannien [9] 
enables the direct evaluation of the crack closure integral (2.2) and 
thus the energy release rate through the use of a finite element 
model. The finite element model starts with the presence of an initial 
crack of length, a. Fig. (2.1a), with tip at node K. The finite element 
solution determines the displacement components, u^, (where u k =>(u k . 
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v k , w k ) of the crack tip node K. An incremental crack extension Aa is 
introduced by replacing the crack tip node K with two separate nodes K 1 
and K'' as shown in Fig. (2.1b). With this new crack geometry taken 
into account, the finite element solution for the nodal displacements 
u k * and uY 1 are found for nodes K' and K' 1 , respectively, under the 
same load. The crack extension is then closed by applying equal and 
opposite forces at nodes K‘ and K ' ' such that their common displacements 
match the displacements found earlier for node K, Fig. (2.1c). These 
nodal forces can be described by 

F k -> < F xk> V F zk> (2 - 3 > 

The energy release rate is then given by [9], 

G(a) = CF xk (u k ' ~ u k ' ' ) +F yk( v k“ v k ’ ' ) +F zk ( w k ‘ “ w k ' ')]/ ?Aa ( 2 - 4 ) 

By resolving the forces and displacements into a "crack coordinate 
system," Fig. (2.2), the respective fracture mode contributions to the 
total energy release rate can be determined. That is, 

Gj (a) = {[F^cosijj-F^si n<}>][cos<J> (w k ' -w^ ' ' )-si mfrfv^ '-v^ 1 ' ) ]} /2Aa (2.5a) 

Gj j (a) = {CFy k c os4>+F z ^si n4>][cos<}> (v^ 1 -v k 1 ' )+si n<t> (w R '-w k 1 ' ) ]} /2Aa (2.5b) 

Sn (a) = t F xk' u k'- u k"X/ 24a 


(2.5c) 
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where Gj, Gjj and Gjjj are the mode-I, mode-II and mode-III contribu- 
tions to the total energy release rate, G, respectively, i.e., 

G (a ) = G j (a ) + Gjj(a) + Gjjj(a) (2.5d) 

2.1.3 Finite Element Considerations in Computing the Energy 
Release Rate 

As previously mentioned, the energy release rate requires the 
evaluation of the nodal forces and displacements necessary to close a 
crack of extended length, a + Aa, back to its original length, a. The 
needed displacements are directly obtained from finite element solutions 
of the initial and extended crack states. Fig. (2.1a) and (2.1b), 
respectively. However, the calculation of the required forces are not 
as obvious. Rybicki and Kanninen [9] computed the forces by placing a 
very stiff "spring" between nodes K 1 and K' 1 , then computed the force 
components in the "spring." This procedure can lead to unnecessary 
approximation errors. An alternative approach will now be presented. 
Consider three separate states. State No. 1 represents the loaded 
initial state. Fig. (2.1a), where node K displaces (u k , v^ ,w k ) . The 
finite element equations (Appendix A) for State No. 1 can be written as, 

- { F x } (2.6a) 

where [Kj] is the global stiffness matrix, {<5^ is the global displace- 
ment vector and { F^ } is the global force vector for State No. 1. Simi- 
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larly. State No. 2, Fig. (2.1b), represents the loaded extended state 
which can be expressed by. 


K 2 ](« 2 ) = {F 2 ( (2.6b) 

and State No. 3, Fig. (2.1c), represents the loaded separated state with 
applied forces to hold nodes K' and K * ' together - which may be written 
as. 


[ K 3 ]{6 3 } = {F 3 } (2.6c) 

Since the forces required to hold nodes K' and K ' ' together are con- 
tained within the {F 3 } vector it is necessary to compute { F^} . Since 
the separated state with applied forces. State 3, Fig. (2.1c), is con- 
strained to displace identical to that of the unseparated state. State 
1, Fig. (2.1a), the displacement vector {<5^} is the same as {fi 3 } with 
the exception of the additional degrees of freedom (i.e., u, v and w 
displacements) for the new node created by separating the crack tip node 
into two nodes. Now, if { <5 ^ } is defined as being the { 6 ^ } vector with 
the additional degrees of freedom it follows that 

{« 3 } = i 6 ! ’ } (2.6d) 

where the additional degrees of freedom are specified as being the same 
as for the initial crack tip node of State 1 since State 3 requires that 
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the displacements of the separated nodes match those of the unseparated 
state. State 1. Note that the new node created by separating the crack 
tip node was numbered as being the crack tip node number plus one and 
all node numbers greater than the crack tip node were re-numbered as 
being one plus the node numbers that they had in State 1. This re- 
numbering procedure guarantees that the half bandwidth will not increase 
by any more than 2 for plane stress or plane strain and by any more than 
3 for generalized plane strain. 

The undeformed mesh of the separated state. State 2, Fig. (2.1b) is 
identical to that of State 3, Fig. (2.1c), and since the stiffness 
matrices do not change for different loading conditions, it follows that 

[k 3 ] = [k 2 ] (2.6e) 

Substituting Eqns. (2.6d) and (2.6e) into (2.6c), the solution 
for {F^} is found to be 


! V l V ! = ,F 3> (2 - 7) 

Hence, for a growing crack problem, the forces necessary to close the 
current crack extension are found by simply multiplying the current 
stiffness matrix by the, modified, previous displacement vector - with- 
out the addition of extra steps or the introduction of unnecessary 
approximation errors. Note that it is not necessary to store the entire 
[k 2 ] stiffness matrix. Since [k 2 ] is a banded matrix, the only contri- 
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butions to the force necessary to close the crack will come from the 
elements containing the crack tip node K. Hence, the elemental contri- 
butions to the force at node K can simply be summed up to give the total 
force at node K. 

2.2 Criteria for Predicting the Direction of Crack Extension 
2.2.1 Modified Griffith Criterion 

The Griffith or energy release rate criterion states that a crack 
will extend when the available energy release rate, G(a), reaches or 
exceeds the critical energy release rate, G c [13]. In a crack problem 
where the crack extension direction is unknown the criterion should be 
modified to state that crack extension will occur in the direction in 
which the available energy release rate, G(a), first reaches the criti- 
cal energy release rate, G c . 

If the critical energy release rate is assumed independent of 
direction then the direction of crack extension can be taken as the 
direction of maximum available energy release rate since this would be 
the direction which would first reach or exceed the critical energy 
release rate. 

Two serious limitations of the modified Griffith criterion, as 
defined above, are that the critical energy release rate may have a 
dependence on the mode of fracture in isotropic problems and it also 
depends on which direction, relative to the material principal coordi- 
nates, the crack extends in fibrous materials. As an example of the 
latter of the two limitations, consider two different mode-I cracks 
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extending in an anisotropic material. The first crack. Fig. (2.3a), 
represents a mode-I crack growing parallel to the fibers and the second 
crack. Fig. (2.3b), represents a mode-I crack growing perpendicular to 
the fibers. Based on surface energy considerations [13], an approximate 
relation for the mode-I critical energy release rate, Gjq, for an iso- 
tropic material is given by 

2 

tto a 

G Tr =— (2.8) 

ic E 


where is the critical applied stress required to cause crack exten- 
sion and E is Young's modulus. Substituting the ultimate strength of a 
composite material, (T300/5208 graphite-epoxy), in the transverse direc- 
tion, Yj, and the modulus in the transverse direction, E£, from Appendix 
C, into Eqn. (2.8), gives an approximate value for G IC for extension 
parallel to the fibers. Fig. (2.3a), that is. 



Similarly, for the crack growing perpendicular to the fibers. Fig. 
(2.3b), 



Hence, the critical energy release rate for the crack growing perpen- 
dicular to the fibers is roughly two orders of magnitude greater than 
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Figure 2.3 Mode-1 Cracks In A Fibrous Material 
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that of the crack growing parallel to the fibers. 

Another disadvantage of the modified Griffith criterion is that it 
requires an additional finite element solution for each possible direc- 
tion that the crack can extend. Referring to Fig. (2.4), for the crack 
defined by oab there are seven possible directions of crack extension, 
from o to c, from o to d, from o to e, from o to f, from o to g, from o 

to h and from o to i. In order to use the modified Griffith criterion, 

seven independent finite element solutions would be required to compute 
the seven possible energy release rates. This is obviously time consum- 
ing and thus a costly procedure. 

2.2.2 Si h Strain Energy Density Criterion 

The strain energy density criterion [2, 15] is based on the local 

value of strain energy density in the vicinity of a crack tip, which is 

direction sensitive. Crack extension is postulated to occur in the 
direction of minimum strain energy density when the strain energy 
density factor, S, (to be defined), attains a critical value, S c . 

For a planar crack in an isotropic material under plane strain. 

Fig. (2.5), the strain energy density in the vicinity of the crack 
tip, -gy-, is given as, [15] 

dU 1/ . 2 _ . . .2 . 2 . 

dV = Trr^ a ll k I +2a 12 k I k II +a 22 k II +a 33 k III^ 

+ non-singular terms 


( 2 . 10 ) 
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where k j , kjj and k j j j are the mode-I, mode-I I and mode-III stress 
intensity factors, respectively. 


a ll = t ^-[( 3 - 4v - cos< ^)( 1+ cos<)))] 
a 12 = ik * 2sin4>[cos<()-(l-2v)] 
a 22 = Tk- 4 ^ 1_v ^ 1-cos ^) + (l +cos ^> )(3 c°siJ>-1)] 

a 33 = 4G 


(2.11a) 


(2.11b) 
(2.11c) 
(2. lid) 


and G is the shear modulus of elasticity and v is Possion's ratio. 

Eqn. (2.10) demonstrates that the strain energy density function 
possesses a (1/r) singularity at the crack tip. Hence the expression 

2 2 2 
S = (3^-^kj+2a^2kjkjj 22^11 33^11 1 ) /"^ 


represents the intensity of the strain energy density field in the 
vicinity of the crack tip. The fundamental hypotheses on crack growth 
in the Sih theory are as follows: 

(1) Crack initiation takes place in a direction determined by the 
stationary value of the strain-energy density factor, i.e.. 



at $ = <{> 0 


(2.13) 


25 


(2) Crack extension occurs when the strain-energy density factor 
reaches a critical value, i.e., 

S c = S(k r k n> k HI ), for <t> = 4> 0 (2.14) 

Exact evaluation of the stresses and strains in the vicinity of the 
crack tip with the current finite element model is not guaranteed since 
there exists a geometric singularity at the crack tip which cannot be 
accurately modeled with the linear-elastic analysis formulated herein. 
Hence, it is not possible to use the Si h theory to determine when the 
crack will extend. However, the Si h theory can be used to determine the 
direction of crack propagation in isotropic materials. 

From continuum mechanics [16] it is possible to write an alterna- 
tive form of the strain energy density at a point in a stressed body, 
i.e.. 


dU 1, . 

jVT = -z[o e +o e +o z +x y +t y +t y ) 
dV 2 V xx xx yy yy zz zz yz yz xz xz xy xy' 


(2.15) 


Neglecting the non-singular terms in Eqn. (2.10) and substituting in the 
expression for the strain energy density factor of Eqn. (2.12) gives. 


dU S 
dV = r 


(2.16) 
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Combining Eqns. (2.15) and (2.16), an alternative equation for the 
strain energy density factor, S, is found to be: 


^ 7^°xx e xx +0 yy e yy + °zz e zz +T yz Y yz +T xz ir xz +T xy Y xy) (2*17) 


There is one serious limitation to the use of the strain energy 
density theory in the current study. This limitation is that the theory 
does not account for the anisotropic strength characteristics of the 
material. Since such properties must be accounted for in fibrous com- 
posites, the Si h strain energy density theory is limited to isotropic 
material applications. It should be noted that the strain energy 
density theory has been used in the past to predict crack growth charac- 
teristics in composite materials [2]. However, the success of such 
studies resulted from assuming that the crack was situated entirely 
within the isotropic matrix between fibers. 

The procedure for implementing the strain energy density criterion 
in the finite element model is briefly described as follows: First, the 

possible directions of crack extension in the model are identified by 
the element sides containing the crack tip node, (node-0 in Fig. 

(2.4)). Second, the stresses and strains are calculated in the adjoin- 
ing elements at the element corners, (points c thru i in Fig. (2.4)). 
Third, Eqn. (2.17) is used to calculate the strain energy density fac- 
tor, S, at the respective points. Last, the crack is assumed to grow in 
the direction in which S is a minimum. 


27 


2.2.3 Tsai-Wu Failure Criterion 

Tsai and Wu [11] postulated that a failure surface in stress space 
exists in the form: 

F.s.+F.jS.Sj =1 i, j = 1.....6 (2.18) 

where F^ and F^- are strength tensors of second and fourth order, 
respectively, and s.. represent a contracted form of the stress tensor 
components in material principal coordinates. For an orthotropic lamina 
under plane stress conditions, Eqn. (2.18) becomes: 


+ X^ s ll + (yT + Y~') s 22~(x~X 
Tc Tc Tc Tc 


ll" v Y t Y ) s 22 + t 2 ^12 


'12 


+2F 12 S 11 S 22 = 1 


(2.19) 


where Xy and Yy represent the tensile strength of the material in the 
fiber and transverse directions, respecti vely , X c and Y c represent the 
compressive strengths, represents the shear strength in the 1-2 
plane and F^ is an interaction term which must be determined from a 
biaxial strength test. 

As a failure theory, the Tsai-Wu criterion has several advanta- 
ges. These advantages include, (1) invariance under rotation of coordi- 
nates, (2) transformation via known tensor transformation laws; and (3) 
symmetry properties akin to those of the stiffnesses and compliances. 
However, for use in this study, it has two serious limitations. The 
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first is that it fails to account for differences in creating new frac- 
ture surfaces at various angles to the material principal coordinates 
and second is that it makes no distinction between tensile and compres- 
sive failure which could cause the Tsai-Wu criterion to choose a com- 
pressive direction of crack propagation. 

The Tsai-Wu criterion was incorporated into the finite element 
model by assuming that the crack would extend in the direction where the 
value of the Tsai-Wu polynomial reached a maximum. Referring to Fig. 

(2.4) , the stresses, o.., are calculated at some fixed distance, r Q , 
away from the crack tip at the various locations dictated by the element 
sides incorporating the crack tip node, (points c thru i in Fig. 

(2.4) ). Next, the values of the Tsai-Wu polynomial were computed at 
these points through the use of Eqn. (2.19), (for the case of orthotro- 
pic plane stress). Last, the crack was assumed to extend in the direc- 
tion for which the Tsai-Wu polynomial reached a maximum. 

The choice of r Q is arbitrary within certain limitations. These 
limitations are that r Q should be greater than zero and less than the 
longest possible path that the crack extension could take and not extend 
through more than one element. Referring to Fig. (2.4), if the possible 
path of crack extension were as shown from node o to c, o to d, o to e, 
o to f , o to g, o to h or o to i , then would be limited to that of 
the segment from o to g since that is the direction of longest possible 
single element crack extension. Note that it would be impractical to 
use a distance greater than og since this would require using stresses 
from an element outside those adjacent to the crack tip. 
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2.2.4 Modified Point Stress and Hashin Criteria 

The criteria for predicting the crack growth direction considered 
up to this point, in their present form, are all unsatisfactory for 
anisotropic materials. As pointed out in Sections 2. 2. 1-3, they all 
fail to account for differences in the work required to create a new 
area of crack surface at different directions in an anisotropic mate- 
rial. Since it is imperative that this distinction be accounted for in 
this study, two new criteria are proposed. The first criterion consid- 
ered is a modification of the point stress criterion of Whitney and 
Nuismer [12], and the second is a modification of the Hashin criterion 
[ 10 ]. 

2. 2. 4.1 Modified Point Stress Criterion 

The point stress criterion of Whitney and Nuismer [12] assumes that 
failure of a notched laminate occurs when the local stress at a certain 
distance, r Q , from the notch tip reaches the strength of the unnotched 
lami nate. 

The modified point stress criterion of this study assumes that a 
crack will grow in the direction of the maximum ratio of normal stress 
to strength at a certain distance, r^, from the tip of an existing 
crack. Note that this is equivalent to assuming that a crack will grow 
perpendicular to the plane of maximum tensile stress in an isotropic 
material. Referring to Fig. (2.6), the normal stress, o ±x ( r,<f>), is 

9 ? 

calculated at some fixed distance, r Q , away from the crack tip. The 

stress, o (r ,4>), is then divided by the tensile strength of the mate- 
99 0 


a** 90-0 



Stress Parameters 
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rial, T iA (<j>), (to be defined), normal to the direction of crack exten- 
ts’ 

sion. This ratio, R (r Q ,<t>) , is then used to predict the direction of 
crack extension by assuming that the crack will extend in the direction 
for which R(r Q ,<}>) reaches a maximum. Quantitatively, the ratio is 
defined as: 


R(r 0 ,$) = 



( r o •♦> 



(♦) 


( 2 . 20 ) 


The value of r Q used in this criterion is subject to the same 
limitations as the value of r Q in the Tsai-Wu criteria. Section 2.2.3. 
That is, r Q should be greater than zero and less than or equal to the 
longest path of possible crack extension while not extending through 
more than one finite element. 


2. 2. 4. 2 Strength, T. ., Along a Given Plane in Anisotropic Materials 

<p<P 

The strength, T (<)>), normal to a given direction, is taken as the 
<p<p 

normal stress required to fail an infinitesimal element of anisotropic 

material along a given plane. In the finite element solution the finite 

element sides dictate the directions of possible crack extension. 

T. was taken in this form to account for differences in the energy 
4><P 

required to create new crack surfaces at arbitrary angles with the 
material principal coordinate system. Such a definition is necessary in 
the proposed model to permit selection of the proper direction of crack 
growth. Further, as shown by Herakovich [7], for example, failure of 
individual lamina in a laminate can occur along planes which are neither 
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parallel to nor perpendicular to the fibers. It appears to be most 
difficult if not impossible to experimentally verify a relation 
for T (9) by testing unidirectional laminates since only one direction 
of crack extension would be present for a given laminate. However, 
three obvious conditions must be met by such an expression. The condi- 
tions are: 

1. For an isotropic material, the strength, T (4> ) , should be constant 

99 

and equal to the ultimate strength of the material, o^, independent 
of 9. 

2. For a crack extending parallel to the fibers in a unidirectional 

composite, T (9) should be equal to the transverse tensile 
99 

strength of the material, Yj. 

3. For a crack extending perpendicular to the fibers, T. ,(9) should be 

99 

equal to the tensile strength of the material in the fiber direc- 
tion, Xj. 

For a unidirectional laminate under plane stress, as shown in Fig. 

(2.7), a simple relationship for 1..U) can be postulated. Removing an 

99 

infinitesimal element at (r Q ,9), Fig. (2.8a), and defining the angle 3 
as the difference between the fiber angle, 0, and the assumed crack 
extension angle, 0, that is 


3=0-9 (2.21) 

It is then possible to isolate yet another infinitesimal element. Fig. 
(2.8b), which gives the orientation of the crack extension relative to 
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Figure 2.7 Unidirectional Laminate With An 
Arbitrary Crack Present 
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the material principal directions. Then, assuming that the strength, 
T^(4>), is balanced by only the transverse and longitudinal strength of 
the material Yy and Xy, respectively, the relation for T (4>) is 
obtained by summing forces in the $ direction as, 

T ( ^ ) = x y sin ^ + Y T cos 2 6 (plane stress) (2.22) 

Testing Eqn. (2.22) against the three conditions. 

Isotropic material: X = Y = a 

T T u 

t H = V 1 "^ + %cos 2 e = o u (2.23a) 

Fracture parallel to fibers: & = 0° 

T = Xysin 2 (0°) + Y T cos 2 (0°) = Yy (2.23b) 

Fracture perpendicular to fibers: & = 90° 

t h = Xysin 2 (90°) + Y T cos 2 (90°) = Xy (2.23c) 

The three conditions specified herein are satisfied. Hence, in princi- 
ple, the expression is acceptable. A plot of T x vs. 4> for plane stress 
is shown in Fig. (2.9) for various values of 6 for T300/5208 graphite- 
epoxy, (properties from Appendix C). The maximums in Fig. (2.9) repre- 
sent the combinations of angles for which a crack grows perpendicular to 
fibers and the minimums represent crack growth parallel to fibers. 
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For a transverse crack extending through a laminate. Fig. (2.10), a 
similar relationship can be derived. Removing the lamina at the crack 
tip and rotating it about the z' axis into material principal coordi- 
nates, Fig. (2.11a), gives the geometry necessary to compute the 

I I I 

strengths T (0), T (0) and T (0) in the x', y', and z' directions, 
xx yy zz 

respectively.. Removing an infinitesimal element, oab, from Fig. 

4 I I 

(2.11a), applying the normal strengths Xy, Yy and T xx ’ in the 1', 2' and 
x' directions, respectively. Fig. (2.11b), and summing forces in the x'- 

I 

direction gives the relation for T xx (0) as. 


T (0) 
xx' ' 


YySi n^B 1 + X-pCOs 2 0‘ 


(2.24a) 


Similarly, removing element, ocd, from Fig. (2.11a), applying 

I l I 

strengths Xy, Yy, and T yy , Fig. (2.11c), then summing forces in the y'- 
di recti on yields 


T yy ( e ) = XySin 2 e + YyCOS 2 0 


(2.24b) 


Since the z' and 3' axis in the lamina coordinate system are the same as 
in the laminate coordinate system. 


T^(e) = Zy (2.24c) 

I 

The strength, 1^(40 » normal to a free edge crack extension. Fig. 
(2.12a), is found by removing an infinitesimal element. Fig. (2.12b), 
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Figure 2.10 Free Edge Crack In 

A Laminated Composite 






applying the respective strengths and summing forces in the direction 
which gives 


= T yy sin 2 4> + t zz cos 2 4> (2.25) 

Substituting Eqns. (2.24a), (2.24b) and (2.24c) into Eqn. (2.25) yields 
that 

7<}>(j, (^ ) = (Xj s i n ^ Q +VjCOS 2 e )sin*> + ZyCOS 2 <f> (2.26) 

Testing Eqn. (2.26) against the three conditions, 

I I I 

Isotropic material: Xy = Yy = Zy = o u 

1 2 2 2 2 
T<M = ( o u s i n e+o u cos 6 )s i n <p + a u cos <j> = a u (2.27a) 

I 

Fracture parallel to fibers: 0 = 0°, <|> = 90° 

T* = (xjsin 2 (0°)+Y|cos 2 (0°))sin 2 (90°) + zjcos 2 (90°) = Y-J- (2.27b) 

I 

Fracture perpendicular to fibers: 0 = 90°, <f> = 90° 

= (XySin 2 (90 0 )+Y T cos 2 (90°))sin 2 (90 0 ) + ZyCOS 2 (90°) = Xy 

(2.27c) 


Again, the three conditions are satisfied so, in principle, the relation 

I 

is acceptable. A plot of vs. <f> for free edge crack extension is 
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shown in Fie,. (2.13) for various values of fiber orientation, 8 , for 
T300/5208 graphite-epoxy , (properties from Appendix C). As in the plane 

I I 

stress case, the values of 4>, (for 0 = 90°), for which T^ is a maximum 

represents crack extension perpendicular to fibers while a minimum 

I 


represents crack extension parallel to fibers. The constant 

I I 

T.. for 8=0° represents matrix mode failure independent of <b. 
<p<p 

I I I 


Last, 


the curve for 8 = 45° has minimums at = Y x , which represents matrix 

4><P l 

I 

mode failure and the maximums never reach Xy since some combination of 
fibers and matrix, thru the width, is always involved for the free edge 
cracks considered. Fig. (2.10). 


2. 2. 4. 3 Modified Hashin Criterion 

The Hashin failure criterion [10] assumes that failure of a trans- 
versely isotropic material will occur in a tensile fiber 
mode, > 0, when: 


(— ) 2 + — (‘uSs) = 1 


(2.28) 


’12 


where Xy is the tensile failure stress in the fiber direction and 1S 
the axial failure shear stress. The Hashin criterion also assumes 
failure to occur in a tensile matrix mode, o ^ + a ^ > 0. when: 


,2 (° 22 + ° 33 ) + _2 ^ a 23 ' a 22 ° 33 ^ + r 2 ^° 12 + ° 13 ^ = 1 ^ 2 * 29 ^ 


23 


12 


43 


r' r 


Or 




Figure 2.13 Free Edge T^vs^For Various Values of Of Q' 
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where Y-j. and are the transverse tensile and shear strength of the 
material, respectively. For the case of plane stress in the 1-2 plane, 
the criterion for fiber mode tensile failure becomes: 


o 2 o 2 
(Y 11 ) + (y 11 ) =1 


(2.30) 


12 


and for matrix mode tensile failure: 


°22 ^ °12 ^ 

("7 / Mr) = 1 
t T ^12 


(2.31) 


where is the shear strength of the material in the 1-2 plane. 


The Hashin criterion does account for distinct differences in fiber and 

matrix mode failure. However, it does not account for an arbitrary 

combination of matrix and fiber mode failure. (In the actual failure of 

composite laminates this feature is necessary since as Herakovich [7], 

for example, has shown failure can occur on a plane which is neither 

parallel to nor perpendicular to the fibers.) 

Hashin [10] proposed that a similar criterion for the failure of 

composite materials could be developed to include the plane on which 

failure occurred, <f> , Fig. (2.14). Such a criterion would predict 

o 

failure when some function g(o.., <J>) satisfied the condition: 

* J 


g(<? i j » ♦) = i 


(2.32) 
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Figure 2.14 Composite Laminate With Modified 

Hashin Parameters 
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and that failure would occur on a plane defined by, 4> q , for which (2.32) 
was first satisfied under monotonically increasing load. 

The modified Hashin criterion of this study assumes that a given 
crack in a composite material will extend in the direction, 4> q , in which 
the left side of Eqn. (2.32) reaches a maximum when evaluated at some 
fixed distance, r , from the crack tip. (Note: r Q is subject to the 

same limitations as in Sections 2.2.3 and 2. 2. 4.1). 

The development of the modified Hashin criterion as used in this 
study is based on developing an expression for Eqn. (2.32). As in the 
case of the point stress criterion, it appears impossible to test such 
an expression experimentally. However, two obvious conditions should be 
met by such a criterion, they are: 

1) For a crack extending parallel to the fibers, the criteria should 
give back the Hashin criterion for tensile matrix mode failure, 

Eqn. (2.31). 

2) For a crack extending perpendicular to the fibers, the criteria 
should give back the Hashin criterion for tensile fiber mode fail- 
ure, Eqn. (2.30). 

Proceeding along the same line as Hashin [10], if the failure criterion 
is taken as. 


+ ( !dL) 2 . !, O > 0 (2.33) 

<J>4> r<j> 

where T , , and T , are the normal and shear strengths and o and o are 
<p<p r<p <?<? r<p 

the normal and shear stresses, respectively, on the plane of crack 
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extension. Fig. (2.14), then the direction of crack extension is given 
by the value of 4> for which the left side of Eqn. (2.33) reaches a 

maximum, provided that a,, is greater than zero. (A o A , less than 

<p<p <M 

zero would represent crack closure.) For the case of plane stress of a 
unidirectional composite laminate, Fig. (2.7), the normal strength, 
T^(<t>), was derived in Section 2.2.4. 2 and is given by Eqn. (2.22). If 
the shear strength, T (<j>), is simply taken as being Sj 2 * 'that is 

T = $12 (plane stress) (2.34) 

Then the failure criterion is complete. Testing the criteria of Eqn. 
(2.33) against the two conditions specified herein, 

(i) Fracture perpendicular to fibers: <j> = 90 + 0, equilibrium of an 

element. Fig. (2.15a) gives 

2 2 

= a yy sin * + °zz cos + " 2T yz sin<t>cos<t> (2.35a) 

2 2 

0 ^ = cos<(>si n<J> (° zz -°yy) + T y Z (cos 4>-si n <fr) (2.35b) 

substituting 4> = 90 + 6, 

2 2 

o H = o^sin (90+8)+o zz cos (90+0 )-2T^ z sin(90+9 )cos(9O+0 ) 

2 2 

= °yy cos 9+0 zz sin 0+2T^ z sin6cos6 


(2.35c) 




and, 
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cos (90+0 )sin (90+6) (° zz°yy^ +T yz [ cos ^ (90+0 ) -si (90+6 ) ] 

-sin0cos6(a Z2 -a^) + t^ z [si n^6-cos^0] (2.35d) 


but, from Fig. (2.15b); 

2 2 

°11 = °yy cos e+0 zz sin 0+ 2 T y Z sin6cos6 (2.35e) 

= -{-si n0cos0 (° zz _a yy) + T yz Csin 2 0-cos 2 0]} (2.35f) 


comparing Eqns. (?.35e) and (2.35f) with (2.35c) and (2.35d), respec- 
tively. 


and 



(2.35g) 



(2.35h) 


Also, for fracture perpendicular to the fibers, Eqn. (2.22) reduces to 
that of Eqn. (2.23c). Substituting Eqns. (2.23c), (2.34), (2.35g) and 
(2.35h) into Eqn. (2.33) yields. 



a 2 

+ (r^ 
*12 


1 


( 2 . 3 5 i ) 


which is precisely the Hashin criterion for tensile fiber mode failure, 
Eqn. (2.30). 
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(ii) Fracture parallel to fibers: <J> = 6 Substituting, <j> = 0 , into 

Eqns. (2.35a) and (2.35b) gives that: 

p p 

a., = a sin 0 + 0 ,_,cos 6-2t sin8cos0 (2.36a) 

<p<p yy zz yz 


and 


2 2 

a p4) = cos0sin0(o zz -a yy ) + x yz (cos 0-sin 0) 

but from Fig. (2.15c): 

2 2 

°22 = °yy sin 0 + °zz cos 6 “ 2T y Z s ' in0cose 


(2.36b) 


(2.36c) 


and 

2 2 

0 j 2 = cos0sin0(o zz -a yy ) + T yz (cos 0-sin 0) (2.36d) 

Comparing Eqns. (2.36c) and (2.36d) with (2.36a) and (2.36b), respec- 
tively, 


and 



(2.36e) 



(2.36f) 


Also, for fracture parallel to the fibers, Eqn. (2.22) reduces to that 
of Eqn. (2.23b). Substituting Eqns. (2.23b), (2.34), (2.36e) and 
(2.36f) into Eqn. (2.33) yields. 



°12 2 
+ i-^) 
12 


1 


(2.36g) 
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which is precisely the Hashin criterion for tensile matrix mode failure, 
Eqn. (2.31). Since the relation developed gives back the Hashin cri- 
teria for tensile fiber mode and tensile matrix mode failure it is in 
principle acceptable. 

2. 2. 4. 4 Finite Element Considerations in Implementing the Modified 
Point Stress and Hashin Criteria 

The steps in the implementation of the modified point stress and 
Hashin criteria are similar to those in implementing the Si h and Tsai-Wu 
criteria. Briefly: 

(i) Determine the elements containing the crack tip node. This gives 

the possible directions of crack extension. 

(ii) Find the minimum element side length. This gives r^. 

(iii) Compute the stresses along the element sides, which define the 

possible directions of crack extension, at r . 

o 

(iv) Use Eqn. (2.20) for the point stress criterion or Eqn. (2.33) for 
the Hashin criterion and assume crack extension in the direction 
which makes (2. 2D), for point stress, or (2.33), for the Hashin 
criterion, a maximum. 


Chapter 3 
RESULTS 


3.1 Isotropic Cases 

3.1.1 Mode-I Crack in an Infinite Plate 

The classical problem of a mode-I crack in an infinite plate. Fig. 
(1.2a), was run as a test of the energy release rate formulation. 
Results for two different finite element meshes, one being much finer 
than the other, were generated for comparison. The computed energy 
release rates were converted to stress intensity factors for ease of 
comparison with theory. 

The boundary condition for the fine mesh, a 306 element x 338 node 
mesh. Fig. (D.l) considered a full crack model assumed specified dis- 
placement loading. Referring to Fig. (3.1a), the boundary conditions 
were: 


at y = 0: v(y = 0,z) = -6 

(3.1a) 

at y = L: v(y = L,z) = <5 

(3.1b) 

at z = 0,2B: traction free 

(3.1c) 


The input parameters; a, 6, L, B, A and B were taken as, 
a = 0.5'', Aa = 0.2a, L = 40a, B = 10a 

A = L/2, 6 = .002a, B = 90° (3.2) 
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The problem was treated as a case of plane strain, constitutive equation 
given by Eqn. (B.18), with material properties listed in Appendix C. 

The crack was assumed to be a virtual crack, i.e., having a width equal 
to zero. 

The coarse mesh, a 68 element x 82 node mesh, Fig. (D.2), assumed 
symmetry about the midplane and also used specified displacement load- 
ing. Referring to Fig. (3.1b), the boundary conditions were: 


at y = 0: 

v(y = 0,z) = -6 

(3.3a) 

at y = L: 

v(y = L,z) = 6 

(3.3b) 

at z = 0: 

Traction free 

(3.3d) 

at z = B: 

w(y, z = B) = 0 

(3.3d) 


The input parameters were taken the same as for the fine mesh, Eqn. 
(3.2). 

The theoretical stress intensity factor, Kj, for a mode-I crack in 
an infinite plate. Fig. (3.1), is given by [17] as, 

Kj = afua) 1 ^ 2 (3.4) 


and the relation between the mode-I stress intensity factor, Kj, and the 
mode-I energy release rate, Gj, for plane strain is given by [17] as. 


K, ■ 


(3.5) 
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where E is Young's modulus and v is Possion’s ratio. Since the theore- 
tical stress intensity factor, Eqn. (3.4), is computed using the remote 
stress, a, the applied displacement, 6, must be converted to an equiva- 
lent stress. From Eqn. (8.18) this relation is given by 

5*«y V = -^ (3 ' 6) 

yy 1-v 

where 



The finite element results for the two meshes considered are pre- 
sented in Table (3.1). The stresses, a, were computed using Eqns. (3.6) 
and (3.7). The comparisons of the predicted stress intensity factors 
with the theoretical stress intensity factors is shown in Table (3.2). 

The results indicated pure mode-I crack extension, i.e., -0 = 0° in 

o 

Fig. (3.1a), as expected. Note that theoretical values for two differ- 
ent initial crack lengths, a, are shown. The first crack length repre- 
sents the actual, initial, crack length and the second represents the 
actual, initial, crack length plus the crack extension. The two crack 
lengths were considered to demonstrate the error incurred in using a 
finite crack extension. That is, the theoretical energy release rate, 
Eqn. (2.1), is based on an infinitesimal crack extension, i.e., limit 
4a * 0, whereas the finite element model introduces a finite crack 
extension, Aa. The results of Table (3.2) indicated that the fine mesh 
gave better results than the coarse mesh and that using the crack length 
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plus the extension in the theoretical stress intensity factor compared 
better with the finite element results. 

3.1.2 Mixed Mode Fracture of an Infinite Isotropic Plate 

The problem of a mixed mode crack in an infinite isotropic plate. 
Fig. (1.2b), was run as a test of the crack growth direction. Results 
were compared with the results predicted by the Si h strain energy 
density theory [2]. The mesh used consisted of 306 elements and 338 
nodes. Fig. (E.l). 

The boundary conditions used were identical to those of the fine 
mesh for the mode-I crack in an infinite plate, i.e., Eqns. (3.1a), 
(3.1b) and (3.1c). Referring to Fig. (3.1a), the input parameters; 
a, 6, L, B, and A were taken as, 

a = "Mr* L = 40a * 8 = 10a 

6 = 0.001", A = L/2 (3.8) 

while the crack inclination angle, 3, was varied from 30 to 90 degrees. 

While simple relations for the mode-I stress intensity factor, Kj, 
and the mode-II stress intensity factor, Kjj, exist for the mixed mode 
problem of Fig. (1.2b), reference [17] points out that there is no known 
relation between the energy release rates and stress intensity factors 
for such a problem. Hence, discussion of the results for the mixed-mode 
crack problem of Fig. (1.2b) will be limited to the crack extension 
di recti on. 


59 


The theoretical crack extension direction, 8 , was computed from 
the Sih strain energy density theory [18], by solving for 6 q from Eqn. 
(3.9), below: 

2(l-v)sin(6 o -20) - {2sin[2(e o -fS)]} - sin26 o = 0 (3.9) 

Finite element results were generated for various crack inclination 

angles, 8, through the use of the modified Griffith criterion. Section 

2.2.1, the Sih strain energy density theory. Section 2.2.2, and the 

modified point stress and Hashin criterion, Section 2.2.4. A plot of 

the theoretical and finite element crack extension direction, 0 , as a 

o 

function of the crack inclination angle, 8, is shown in Fig. (3.2). 

The results for the modified Griffith criterion were generated by 
evaluating the value of the crack closure integral, through the use of 
Eqn. (2.4), for all possible paths of crack extension present, e.g., see 
Fig. (2.4), and assuming that crack extension would occur .in the direc- 
tion of a maximum energy release rate. The results consistently predic- 
ted crack extension in a direction in which the mode-I energy release 
rate, Gj(a), made up 99% or more of the total energy release rate. 

The results for the Sih strain energy density criterion were evalu- 
ated using Eqns. (2.13) and (2.17). The strain energy density factor 
was considered as the sum of two components [18], one due to a change in 
volume, Sy, and one due to a change in shape, S d , 


(3.10) 


-do (degrees) 
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Figure 3.2 Crack Extension Angle,-^, vs. Crack 
Inclination Angle, (3 , For Isotropic Mixed Mode Fracture 


The results always predicted crack extension in a direction in 
which S v > > S^. This conforms with the concept [18] that fracture 
occurs along a plane where S y > S^. 

The modified point stress cases were evaluated using Eqns. (2.20) 
and (2.22) with 


X T = Y t = a u (3.11) 

The crack was assumed to extend in the direction for which, R(<f>), Eqn. 
(2.20) was a maximum. 

The results for the modified Hashin criterion were generated 
through the use of Eqn. (2.33), with 



(3.12) 


Crack extension was assumed to occur in the direction which maximized 
the left hand side of Eqn. (2.33). 

All of the tested theories gave identical results, as indicated in 
Fig. (3.2). The small differences between the theoretical and finite 
element predicted values was attributed to there being only a finite 
number of crack extension paths available in the finite element model 
compared to an infinite number in the analytical Sih theory. However, 
the finite element model always predicted crack extension along the 
closest available crack extension path to that of the theoretical crack 


extension direction. 
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3.2 Off-Axis, Unidirectional Composite Tension Specimens. 

A tension test of a unidirectional coupon with a small crack pre- 
sent, Fig. (1.3a), was simulated as a test of the crack extension direc- 
tion. The mesh used consisted of 306 elements and 338 node points. Fig. 
(D.l) . Referring to Fig. (3.1a) and (2.7), the boundary conditions were 
chosen to simulate the grips of a tension test machine and are given by: 

at y = 0: v(y = 0,z) = -6, w(y = 0,z) = 0 (3.13a) 

at y = 2L: v(y = L,z) = <5, w(y = L,z) = 0 (3.13b) 

at z = 0,2B: traction free (3.13c) 

The dimensions were chosen to be similar to a typical tensile coupon, 
(Reference [19] suggests L > 308 be used.), and were taken as, 

B = 0.25' ' 

L/2B = 15 

The applied displacement load was chosen to be, 

6 = 0.001" = 3750 (in) (3.14c) 

The initial half-crack length was chosen as. 


(3.14a) 

(3.14b) 


a 


0.005 in 
sinB 


(3.14d) 
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and the other parameters, the fiber direction, 0, the crack inclination 
angle, 8, and the crack position. A, were varied to simulate various 
conditions. The problem was treated as being a case of orthotropic 
plane stress with the constitutive relation being that of Eqn. (B.13). 
The material system chosen was T300/5208 graphite-epoxy, properties 
listed in Appendix C. The crack was assumed to be a virtual crack, 
i.e., having zero initial width. 

3.2.1 Comparison of the Crack-Extension Direction Theories for a 30° 
Lamina 

The theories for predicting the crack extension direction, i.e., 
the modified Griffith theory. Section (2.2.1), the Si h strain energy 
density theory. Section (2.2.2), the Tsai-Wu theory, Section (2.2.3), 
the modified point stress theory. Section (2.2.4), and the modified 
Hashin theory. Section (2.2.4), were compared against the expected crack 
extension path for a 30°, off-axis, unidirectional tensile specimen with 
a crack orientated along the z axis, i.e., 8 = 90° in Fig. (3.1a). 

The predicted crack extension path for the modified Griffith theory 
is shown in Fig. (3.3a), the Si h strain energy density criterion in Fig. 
(3.3b), the Tsai-Wu theory in Fig. (3.3c), the modified point stress 
theory in Fig. (3.3d) and the modified Hashin theory in Fig. (3.3e). 

The experimentally observed direction is shown in Fig. (3.3f). Experi- 
mental results for graphite-epo*y [6, 7 and 8], Boron-aluminum [20], and 
graphite polyimide [21], unidirectional composites all indicated that 
failure of unnotched specimens and fracture of notched specimens 
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occurred along planes parallel to the fibers. Fig. (3.3a) thru (3.3e) 
indicate that the modified point stress and the modified Hashin theories 
were the only theories that predicted the correct crack extension path. 

The reason that the modified Griffith, the strain energy density 
and the Tsai-Wu theories predicted incorrect crack extension paths was 
because none of these theories account for the differences in the energy 
required to create crack extension surfaces at arbitrary angles to the 
fi bers. 

Since the modified point stress and Hashin criteria were the only 
criteria to yield accurate results, further case studies were limited to 
these two criteria. 

3.3.2 Variation of the Modified Point Stress and Hashin Functions for a 
30° Lamina. 

The modified point stress function, R ( r Q ,4> ) , was given by Eqn. 
(2.20) as 


where a was taken as the normal stress and T xx , defined by Eqn. 

<P<P <p<p 

(2.22), the strength normal to a plane of possible crack extension. The 
modified Hashin function can be defined as H(t q ,<{>), where from Eqn. 
(2.33), 

2 a . 2 

H(r 0 ,*) = (-**■) + H 4 ) 

T H ^"r4> 


(3.16) 
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where was the shear stress and T the shear strength, Eqn. (2.34), 

on the plane of possible crack extension. 

Finite element predictions of R(r Q ,4>) and H ( r o ,4> ) were computed as 

a function of 4> for three different values of r for the 30° lamina of 

o 

Section 3.2.1. (The finite element mesh, the boundary conditions and 

the crack geometry were identical to those of the problem considered in 

Section 3.2.1). The values for the R ( r Q ,4> ) function were normalized to 

their maximum values and plotted in Fig. (3.4) for three different 

normalized values of r . The r values were normalized with respect to 

oo 

r max , where r max was the limiting value of r for the finite element 
o’o o 

mesh used, (e.g., see Section 2. 2. 4. 4). The values for the H(r Q ,<J>) 
function of the modified Hashin theory were also normalized to their 
maximum value and were plotted in Fig. (3.5) for the same three norma- 
lized r values, 
o 

The results of Fig. (3.4) and (3.5) achieve maximum values at 4> = 

210° for all of the various r /r max values. Since the criteria both 

o o 

assume crack extension in the direction where R(r Q ,4>) or H ( r 0 ,4> ) reach a 

maximum, both criteria choose the expected value of 4> = 210°, (e.g., see 

Section 3.2.1). Figure (3.4) and (3.5) also indicate that the trends 

remain unchanged though the values differ slightly for the three r 0 /r™ ax 

values. This indicates that the prediction of the crack extension 

direction is fairly insensitive to the value chosen for r . Note that 

o 

in the case of the modified point stress theory. Fig. (3.4), the maximum 
value of R/R max increased with increasing distance from the crack tip. 
This trend was just the opposite of what was expected since the stresses 
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Figure 3.4 Normalized Modified Point Stress Ratio vs.0 
For A 30 ? Unidirectional, Laminate 
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0 (degrees) 

Figure 3.5 Normalized Modified Hashin Function vs.0 
For A 30? Unidirectional , Laminate 
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are actually singular at the crack tip, i.e., at r Q = 0. The reason for 

this discrepency is that as r Q /r™ ax approaches unity, the stresses are 

computed nearer the element corners which introduces numerical errors in 

the stress computation. This stands as a good example of why such 

criteria are limited to the prediction of crack extension direction and 

not the load which would cause crack extension. The values for the 

normalized modified point stress ratio. Fig. (3.4) and the normalized 

modified Hashin ratio. Fig. (3.5), for 0° <<{><= 180° were not shown 

because of numerical difficulties in the computation of the stresses in 

this regime. However, the modified point stress results seemed to 

indicate a slight increase in ^ ( r Q ,4> ) for 0° < <{> < 30° and some negative 

values between 4> = 30° and $ = 180°. Similarly for the modified Hashin 

ratio. Fig. (3.5), the values forO° < <j> < 180° seemed to indicate slight 

increases in H (r^ ,4> ) around <t> = 30° and <j> = 120° with some negative 

values for o between <f> = 30° and = 180°. Comparison of the modified 
<P9 

point stress ratio, R( r »$), Eqn. (3.15), and the modified Hashin 
ratio, H (r , 4>) , Eqn. (3.16), reveal that 

H(r ,♦) = [R(r 0 ,4>)] Z + (y^) (3.17) 

r<J> 

Hence, the differences in Fig. (3.4) and (3.5) represent the effects of 
squaring the point stress ratio, R ( r Q ,4* ) , plus the shear effect of 
including the square of the shear ratio, a /T . 
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3.2.3 Crack Extension in Unidirectional Laminates 

The modified point stress and the modified Hashin criteria were 
used to predict the crack extension direction and the modified crack 
closure technique was used to predict the energy release rate for 
several different unidirectional laminates. The reason the study of 
crack extension in unidirectional laminates was limited to the modified 
point stress and Hashin criteria was because they were the only criteria 
which accounted for differences in the energy required to create new 
fracture surfaces at arbitrary angles to the material principal system — 
as demonstrated for a 30° lamina in section 3.2.1. The mesh, boundary 
conditions and material system were identical to those explained at the 
beginning of section 3.2. A virtual crack, i.e., one having zero ini- 
tial width, was assumed. The orientation of the crack was assumed to 
be 0 = 90° and A = L/2, Fig. (3.1a). 

The results are shown in Table (3.3) where 6 is the material prin- 
ciple coordinate system orientation, <)> o the orientation of the plane of 
crack extension, G is the total energy release rate, Eqn. (2.4), and the 
% mode-1 and % mode-II values represent the mode-I and mode-II contri- 
butions to the total energy release rate, G. The expected crack exten- 
sion directions were based on experimental results for graphite-epoxy 
[6, 7], Boron-aluminum [20], and graphite polyimide [21] composites. 

The results of these experiments all indicated that failure of unnotched 
specimens and fracture of notched specimens occurred along planes par- 
allel to the fibers and that <j> was always greater than 180° except in 
the case of 9 = 0°, where extension could occur in either a <(» = 0° 


Table 3.3 Finite Element Crack Extension Predictions 
For Unidirectional Laminates 
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or 4> o = 180° direction. The results of Table (3.3) indicate that both 
the modified point stress and Hashin criteria choose the correct path of 
crack extension in all of the cases considered. The fact that both 
criteria also choose dual values of <t> = 0°, 180° for the 0 = 0° speci- 

mens was also promising since this is what would be expected. While no 
direct comparison of the computed energy release rates could be made, 
the computed values were similar to and in the same range as experimen- 
tal values obtained by Wang and Crossman [6] for double side notched 
graphite-epoxy specimens. The % mode-I and the % mode-II values were 
also compared with the results of Wang and Crossman [6] and were found 
to exhibit the same trends and range of results. 

3.2.4 Effects of Crack Orientation on the Fracture Characteristics of 
Unidirectional Laminates 

The finite element model was used to predict the fracture charac- 
teristics of a 30°, unidirectional laminate. The cases considered were: 

1. Influence of crack position, referring to Fig. (3.1a), a crack near 
the grip of the tensile machine was simulated by specifying 

that £ = 0.1 with B = 90°. 

2. Influence of specimen aspect ratio, referring to Fig. (3.1 a), a 
short specimen was modeled having length L/2B = 2 with B = 90°. 

3. Influence of crack orientation, referring to Fig. (3.1a), specimens 
with B = 30°, B = 60°, and B = 90° were considered. 
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All of the cases assumed that a virtual crack existed and used the same 
mesh, boundary conditions, loading, material and geometry, unless stated 
otherwise above, as described at the beginning of section 3.2. 

The results of case 1, Table (3.4), indicated that the crack direc- 
tion for the crack near the grip, A/L = 0.1, extended in the same direc- 
tion, <(> , as did the center cracked specimen, A/L = 0.5, while the 
o 

energy release rates, G, and the % mode-I and % mode-II contributions to 

the total energy release rate differed by 10-20%. 

The case 2 results. Table (3.5), indicated that the crack extension 

direction, <t> , was insensitive to the specimen aspect ratio and that the 
o 

resulting energy release rates differed by about 20%. However, the type 
of fracture that occurred was the opposite for the two cases. That is, 
the long specimen, L/2B = 15, fractured in a mainly mode-I fashion while 
the short specimen, L/2B = 2, fractured in a mainly mode-II fashion. 

This trend is not surprising since Nemeth, et al. [19] have shown that a 
completely different stress state exists in short specimens, i.e., L/2B 
= 5, as compared to long specimens, i.e., L/2B = 15. 

The results of case 3, Table (3.6) predicted the same crack exten- 
sion directions, , and similar energy release rates and fracture modes 
for the three crack inclination angles, 6, considered. These results 
are consistent with the experimental results of [6], [7], [20], and [21] 
which indicate that fracture of unidirectional laminates always occurs 
on planes parallel to the fibers. 


Table 3.4 Influence Of Crack Position On The Fracture 

Characteristics Of 3 0? Unidirectional , Laminates 
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Table 3.6 Influence Of Crack Orientation On The Fracture 
Characteristics Of 30 °, Unid motional, Laminates 
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3.3 Free Edge Crack Growth in Laminated Composite Tensile Specimens 
Cracks located at the free edge of laminated tension specimens. 

Fig. (2.10), were considered. The crack extension paths and the result- 
ing energy release rate vs. crack length curves were considered. 

The solution procedure consisted of obtaining a generalized plane 
strain solution of the front face, i.e., the y-z plane, under an applied 

load of e then, using the results of the front face solution as 
xx 

applied loads, obtain a fracture mechanics solution to the free edge, 
i.e., the y'-z 1 plane, via subsequent generalized plane strain solu- 
tions. (The generalized plane strain formulation is described in Appen- 
dix A.) 

The generalized plane strain solution of the front face was 
obtained using a 132 element by 150 node finite element mesh. Fig. 

(D.3). Referring to Fig. (2.10) and Eqn. (A. 4), the boundary conditions 
and geometry used were: 


B • 

= 0. 

25' ' 

, H = 0.02", £ = 0.001 

(3.18a) 




XX 

at 

y = 

0 , 

U(y = 0,z) = V(y = 0,z) = 0 

(3.18b) 

at 

y = 

B, 

traction free 

(3.18c) 

at 

z = 

0 , 

W(y,z = 0) = 0 

(3.1 8d ) 

at 

z = 

H, 

traction free 

( 3 . 18e ) 


where quarter symmetry was assumed. 


Since the generalized plane strain solution assumes stresses and 
strains to be independent of the out of plane coordinate, the study of 
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the free edge crack growth problem can be modeled exactly only when 

there is no variation of these quantities in both the x and x* direc- 

tions. Since the x-direction is the out of plane direction. Fig. 

(3.6a), for the front face problem and the x'-direction is the out of 
plane direction, Fig. (3.6b), for the free edge problem, the study is 
limited to cases which exhibit x and x‘ independence. The only lami- 
nates which exhibit this quality are uni di rectional laminates because of 
the absence of edge effects. However, if the thru the thickness varia- 
bles of the free edge problem. Fig. (3.6b), are assumed constant and 
equal to the values at the free edge of the front face problem, i.e., at 

y = B in Fig. (3.6a), then the corresponding free edge boundary condi- 

tions become Fig. (2.10), (3.6a), (3.6b), 


at y' = 0, W (y ' = 0,z') = W(y = B,z) (3.19a) 

V ' (y ' = 0,z ' ) = -[U(y = B,z) + e xx -(x = C)] (3.19b) 

at y' = b, W (y ' = 2C,z') = W(y = B,z) (3.19c) 

v (y * = 2C , z 1 ) = -[U(y = B,z) + c^(x = -C)] (3.19d) 


where 2C is the free edge modeled length. 

at z' = 0, W' (y ' ,z ' = 0) = 0 (3.19e) 

at z' = H, traction free ( 3 . 19 f ) 
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with: 


u'(y'.z') = o 

(3.19g) 

xJcfz'l - T»(y - B >*) 

(3 . 19h ) 

0.' =0. - 90° 
1 1 

(3.191 ) 


Several laminate configurations were considered as candidates for the 
analysis of free edge crack growth. Tests were preformed to find which 
laminate configurations could best be modeled as generalized plane 
strain problems on both the front face and the free edge. The tests 
consisted of first obtaining the front face solution then applying the 
corresponding boundary conditions, Eqns. (3.19a-i), to an uncracked free 
edge model. The results of the tests were analyzed to determine which 
laminates gave the best correspondence between the stresses and strains 
at the free edge as well as approximated independence of the displace- 
ments in the x and x' directions, Fig. (2.10), for the front face and 
free edge problems. The results indicated that both angle-ply and 
cross-ply laminate configurations gave a reasonable correspondence of 
free edge stresses and strains. However, only the cross-ply config- 
urations approximated independence of both x and x'. Fig. (2.10), under 
the applied load. Hence, the analysis of free edge crack growth in this 
study was limited to cross-ply laminates. 

Experimental results [6, 22] indicate that laminates which contain 
90° plys along with other plys where 6 ^ 90° can exhibit transverse 


81 


crack growth in the 90° plys and delaminations between the 90° and 6 j- 
90° plys. As an example of the quasi-three dimensional capabilities of 
the finite element model the nature of crack growth in a [QO^/O^ls, 
T300/5208, graphite-epoxy tensile specimen was considered. Referring to 
Fig. (3.7), both transverse crack growth and delamination crack growth 
were considered and the results of the two analysis were compared. The 

method of analysis consisted of testing two cases. Case 1 was to trace 

the crack growth of an initial transverse crack and case 2 was to trace 

the crack growth of an initial delamination. Both cases used the boun- 

dary conditions of Eqns. (3.19a-c). The transverse crack case used the 
crack geometry of Fig. (3.1b) and the delamination case used the crack 
geometry of Fig. (3.8). Both cases were modeled as the free edge of an 
8-ply tensile specimen with 

B = 0.25in, H = O.Olin, C = 2H, a = 0.77H (3.20) 

and an applied normal strain of. Fig. (2.10), 

e = 0.001 (3.21) 

The modified point stress criterion. Section 2. 2. 4.1, was used to pre- 
dict the direction of crack extension and the modified crack closure 
method was used to compute the energy release rates. 

The resulting crack growth sequence for the initial transverse 
crack, case 1, is presented in Fig. (3.9). The case 1 results indicated 
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that the transverse crack extended through the 90° plys down to the 0° 
plys when it turned into a delamination. Note that this does not neces- 
sarily indicate that a transverse crack will turn into a delamination 
but just that, if the transverse crack was the only crack present and if 
it was to continue to extend that it would extend as a delamination 
before extending through the 6=0° plys. The plot of the corresponding 
energy release rate vs. crack length is shown in Fig. (3.10). Point A 
represents the initial crack length and point B represents the point at 
which the transverse crack turned into a delamination. The critical 
energy release rate for such a crack, is given by Wang and Crossman [6] 
for a mode-I transverse crack in graphite-epoxy as 

G IC » 0 . 9 — -- ip -- (3.22) 

in 

Since the available energy release rate curve. Fig. (3.10), was less 

than 0.6 ^ n ’v,— s for the entire crack extension sequence, the crack would 
in 

not have extended at the applied load level of e xx = 0.001. However, if 
the load were increased the G vs. a curve would have translated up, 
retaining the exact same shape, [6], C-D in Fig. (3.10) for example, and 
would have predicted crack growth until the crack length was such that 
the available energy release fell below the critical energy release 
rate. Regardless of applied load, however, the shape of the G vs. a 
curve. Fig. (3.10), indicates that G is decreasing with increasing a. 
Hence, there exists a certain range of e values for which the trans- 
verse crack would grow then arrest when G falls below G c . This trend is 



Figure 3.10 G vs. a For Transverse Crack Growth 

In A Cross-Ply Laminate 
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supported by many experimental results of [6, 22], for example, which 
show that transverse cracks will grow through the 90° plys, arrest and 
form other transverse cracks at regular spacings. 

The results for the initial delamination crack, case 2, indicated 
that the initial delamination would immediately turn into a transverse 
mode. In comparing the point stress ratios of the initial delamination 
crack with the initial transverse crack at the point it turned into a 
delamination, it was found that the point stress ratio was roughly two 
times larger for the delamination turning into a transverse crack than 
for the transverse crack turning into a delamination. This indicates 
that if both a delamination and a transverse crack were present in a 
given specimen, the delamination would turn and grow into a transverse 
crack before the initial transverse crack would extend. This trend is 
supported by the experimental results of Harris and Orringer [22], which 
indicate that transverse cracks can branch off from delaminations. 


Chapter 4 
CONCLUSIONS 

The present investigation has been concerned with predicting the 
direction of crack extension as well as the load to cause extension in 
composite materials. The results of the present study indicate that 
failure criteria can be used to predict the direction of crack extension 
and that an energy release rate approach, implemented through the use of 
a modified crack closure integral, can be used to determine when a crack 
extends and if crack arrest will occur. The finite element models 
presented herein were formulated for two-dimensional and quasi three- 
dimensional analysis. However, the procedures and methods developed can 
be applied to full three-dimensional analyses as well. 

It was found that criteria for predicting the direction of crack 
extension should account for differences in the energy required to 
create crack surfaces at arbitrary angles to the material principal 
system. The Griffith criterion, the Tsai-Wu failure criterion and the 
Si h strain energy density theory all were unsatisfactory in this regard; 
but the modified point stress and Hashin criteria provided good predic- 
tions for crack growth direction. It was also found that since the 
procedure developed herein assumes that crack extension will occur along 
the element sides adjacent to the crack tip node, an incorrect direction 
of extension can be chosen if no element sides coincide with the actual 
direction of extension. However, the model always chooses the closest 
direction available to the actual direction of extension. This problem 
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can easily be overcome by first finding an approximate direction of 
extension using a coarse mesh then refining the mesh in the area of the 
direction chosen by the coarse mesh. 

The present investigation has also shown that future research is 
warranted in the following areas: 

1. The extension of the current work into a full three-dimensional 
model. A full three-dimensional model should be formulated to 
account for the three-dimensional crack growth characteristics of 
many laminates. 

2. Experimental work in the area of critical energy release rates. 
Experimental work should be performed to determine the effect on the 
critical energy release rate of cracks extending at arbitrary angles 
to the material principal system so that the load to cause failure 
can be accurately defined. 
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APPENDIX A 


Linear Elastic Finite Element Relationships 

The finite element used is shown in Fig. (A.l). The element is a 
four node, general quadralateral , isoparametric element. The element 
uses linear interpolation as described by Segerlind [23]. The details 
of the finite element concept are also given in reference [23], The 
technique involves mapping a distorted shape in the Cartesian (y,z) 
Coordinate System into a square in the Local (S,n) Coordinate System 
where £ and n range from -1 to +1. The relationship between the global 
Cartesian and the local coordinates is 

Y = N 2 Yj + N 2 Y 2 + N 3 Y 3 + N 4 Y 4 = i = 1,4 

Z = NjZj + N 2 Z 2 + N 3 Z 3 + N 4 Z 4 = i = 1,4 (A.l) 

where the N. (£, n) are the interpolation functions for the four node 
points and Y^ , Z-j are the Cartesian coordinates of the nodes. The inter- 
polation functions are given by 

Nj = }(l-€)(l-n), N 2 = {(lH)(l-n) 

N 3 = }(l+5)(l+n), N 4 = |(l-C)(l+n) (A.2) 

For an isoparametric element, the same interpolation functions are used 
for the assumed displacements as for the geometry. Hence, for plane 
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Figure A.1 Global And Local Coordinate System - 
4 Node , Isoparametric Element'' 
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stress or plane strain, where there exists two unknown degrees of free- 
dom per node, the interpolation becomes 

v(y,z) = N i v i i = 1,4 

w(y,z) = N i w i i = 1,4 (A. 3) 

For generalized plane strain, where there are three degrees of freedom 

per node, the interpolation is given by 

u(x,y,z) = U(y,z) + e xx *x = N i u i + e^x i = 1,4 

v(x,y,z) = V(y,z) = N i v i i = 1,4 (A. 4) 

w(x,y,z) = W(y,z) = N 1 m 1 i = 1,4 

where u, v and w are the x, y and z displacements, respectively. u n - , v^ 
and w-j are the unknown values at the ith node and are functions of y and 
z only, and e is the total strain in the x direction, which is assumed 
to be constant and either known or unknown. 

The strain-displacement relationships are derived based on the 
small strain - small displacement theory. For the three dimensional 
(generalized plane strain) case, these relationships may be written as 



And for plane strain or plane stress these relationships become 



Substituting (A. 4) into (A. 5), for generalized plane strain. 



where 

I is the 3x3 identity matrix 

{q} is the 12 x 1 vector of nodal displacements 

given by 
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{q} 


r "\ 
U 1 

V 1 

W 1 

u 2 


v 2 

w 2 

u 3 

v 3 

w 3 


> 


u 4 


v 4 


IN 
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(A. 8) 


Substituting (A. 3) into (A. 6), for plane stress or plane strain. 


where 




3 

e 

yy 


ay 0 

e 

_ 

o 5- 

' zz 


9z 

3 3 

V 


iz "sy 





CIN 1 |IN 2 |lN 3 |IN 4 ]{q} 


(A. 9) 


I is the 2x2 identity matrix 
{q} is the 8x1 vector of nodal displacements given by 
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{q} - 


'V 

W 1 

v 2 

w 2 

v 3 

w 3 

v 4 

_ w 4. 
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(A. 10) 


The [B] matrix (strain-displacement relationships) for generalized plane 
strain can be defined as 


xx 

0 

0 

{e} = [B]{q} + 0 

(6x1) (6x12) (12x1) 0 

0 


(A. 11) 


where by comparing (A. 11) and (A. 7), for generalized plane strain. 
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[B] = 


3 

17 

0 

0 

0 

3 

17 

3_ 

9 y 


0 

3 

■57 

0 

3 

■57 

0 


3 

3x 


0 

3 

■57 

3 

5y 

3 

17 


— — n 


[INj | IN 2 I IN 3 | IN 4 ] 


(A. 12) 


For plane strain or plane stress, the [B] matrix can be written as 


{«} = CB]{q} 
(3x1) (3x8) (8x1) 


(A. 13) 


Comparing (A. 13) and (A. 9), for plane stress or plane strain. 


[B] = 


37 

0 

0 

3 

17 

3 

3 

9z 

■57. 


[IN 1 |IN 2 |IN 3 |IN 4 3 


(A. 14) 


Recall from (A. 2) that the are functions of the local coordinates 
K and n. In order to determine the elements of the [B] matrix a rela- 
tionship between the derivatives in the global (y,z) and the local 
(5, n) coordinate systems is needed. This relationship is given by the 
Jacobian matrix [J] of the transformation where 


100 


3N. 

1 

"Sy 

3N. 

1 

■57 


= [J] 


-1 


and 


3N 


35 


<1 


3N. 

T5y 


[J] = 


CJ] 



f 3 y 

■ST 

3z* 


h. 

L3n 

3 Z 

3n_ 

(A.l) 

into 

(A. 15) 

’ 3N 1 

35 

3N 2 

35 

3N 3 

3Ni 

w 

3N 2 

TrT 

3n 
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3N < 

35 

3N, 


3n 


(i = 1, 2, 3, 4) 


*1 

y2 

*3 

y 4 


2 1 

z 2 

z 3 

z 4 


(A. 15a) 


(A. 15b) 


(A. 16) 


The stress-strain relationships are derived in Appendix B. For 
general purposes they can be written as 


M - CD]{e} 


(A. 17) 


For generalized plane strain, {a} and { e } are 6x1 vectors, and are 
given in Appendix B by Eqns. (B.9) and (B.10), respectively. 

For generalized plane strain, the [D] matrix takes on the values of 
the [C] matrix as given by Eqn. (B.8). For plane stress or plane 
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strain {a} and {e} are the 3x1 vectors given by Eqns. (B.14) and 
(B.15), respectively. For orthotropic plane stress the [D] matrix takes 
on the value of the 3x3 [0] matrix of Eqn. (B.13), while for isotropic 
plane strain the [D] matrix becomes the 3x3 [0] matrix of Eqn. (B.18). 

The total potential energy, n, of a given finite element is the sum 
of the strain energy, U, and the work of external loads, W. The strain 
energy, U, of the element is 

U = i / {o }{e } T d V (A. 18) 

V 

and the work of external loads, W, is 

W = -{q)lF} T (A.19) 

where {F} is appl ied mechanical load, (traction), vector. 

Hence, the total potential energy of the element is given by 


n = U + W 


(A. 20) 


Substitution into (A. 20) in terms of the matricies and vectors described 
herein and minimizing with respect to the unknowns yields the finite 
element equations 


[K]{ql = {f| 


(A. 21) 
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where [K] is the elemental stiffness matrix given by 

DC] = / [B] T [0][B]dV (A. 22) 

V 

when mapped into the S,n coordinate system 

dV = |J|d£dn (A. 23) 

where |J| is the determinate of the Jacobian matrix. The limits of 
integration are -1 to +1 in both £ and n. Hence, Eqn. (A. 22) becomes 

1 1 T 11 

[K] = / J [B] [D][B] | J | d£dn = J J G(£,n)d£dn (A. 24) 

-1 -1 -1 -1 

In order to evaluate the stiffness matrix, [K], a numerical integration 
is necessary. Using a 2 x 2 Gauss rule Eqn. (A. 24) can be evaluated as 

2 2 


[K] = Z Z HHG( a b.) 
i =1 j=l 1 J 1 J 

(A. 25) 

where G(£,n) = [B] T [D][B] j J | 

(A. 26) 


(a.j,bj) are the coordinates of the four Gauss points given by 



(A. 27) 
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and (H^.Hj) are the corresponding weight functions given by 


Hi = 1,1 

H j " 1.1 


(A. 28) 


The elemental matricies are then assembled into a global system of 
equations, the prescribed boundary conditions are imposed and the system 
of equations is solved for the unknown displacements. (This procedure 
is explained in most finite element books [23].) 

The strains can be found at any 5 ,n location within an element 
through the use of Eqn. (A. 7) for generalized plane strain and through 
Eqn. (A. 9) for plane stress or plane strain. The strains can be conver- 
ted to stresses by using the stress strain relation of Eqn. (A. 17) where 
[D] takes on the values of the appropriate constitutive relation (Appen- 
dix B). 
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APPENDIX B 

CONSTITUTIVE RELATIONS 

The constitutive relation for an orthotropic material in the prin- 
cipal material directions, Fig. (B.l), is given by Jones [24] as 

= [C]{e}j (B.l) 

where 
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(B.4) 
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V +V V 

12 32 13 
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'22 33 
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V 13 +V 12 V 23 

13 = E n E 22 4 


1 — V V V +V V 

r - 13 31 23 21 13 r 

^22 E E A * l 23 " E F A * ^ 
11 33 t ll t 22 


1 - V 12 V 21 


(3.5) 


c 44 - g 23' c 55 “ g 13- c 66 m G 12 


-V ^ V -V V -V V -?V V V 

12 21 23 32 31 13 21 32 13 

E 11 E 22 E 33 


and 







(B.6) 
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For a 9 rotation about the 3, (z), axis. Fig. (8.1), the consti- 
tutive relation becomes 


where 


[C] - 


M = CC]{e} 





(B.7) 


(B.8) 


(B.9) 


(B.10) 
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If m = cos 0 and n = sin 0 then the components of the [C] matrix are 
given by 


'11 


'12 


'13 


'16 


'22 


'23 


'26 


4 2 2 4 

m + 2m n (Cj 2 +2Cgg) + n C 22 

m n (Cj j + C 2 2 — 6 6 ^ + ^ +l "* ^12 
2 2 

m C i 3 + n C 23 

mn [m (C^ j 12”^^66 ^ + *”* ^12""^'22 +2 ^'66^ 

n 4 C n + 2m 2 n 2 (C 12 +2C g6 ) + m 4 C 22 

2 2 

n C 13 + m C 23 

mn[n (C ^ j ""^12 _2 ^66 ^12”^22^ 2 ^66^ 


'33 


= C 


33 


C 36 = mn ^ C 13~ C 23 ^ 


'44 


m C 44 + n C 55 


£45 = mn ( c 55- c 44) 


'55 “ n C 44 + m C 55 

'66 = m n 2?) + ^66^ m _n ) 


(B.ll) 


For an orthotropic material under a state of plane stress with 
a 0 rotation about the 3, (x) axis, Fig. (B.2), Jones [24] gives the 
constitutive relation as 


{ 0 } = CQ] { e } 


( B - 1 2 ) 
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Figure B.2 Material And Global Coordinate System 
For A Two-Dimensional Lamina 


no 


where 
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0 i 2 
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®12 
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( B .13) 


(B.14) 


(B.15) 


and the Q.. terms of Eqn. (B.13) are given in terms of the C. . terms of 
1 J 1 J 

Eqn. (B.ll) by 




C i 3 C j 3 



(B.16) 


For an isotropic material under plane strain, Frederick and Chang 
[16] lists the constitutive relation as 


f«l - [c]H 


(B.17) 


Ill 


ORiGIS'tAi. r --vVr *_l i -■ 
OF POOR QUALITY 


where 


[ 5 ] 


E(l-v) 

(l+v)(l-7v7 




0 


1 0 


TTr-vy 


(B.18) 


W} 






(B.19) 


(B .20) 


and 


a = v(ct +o ) 
zz v yy zz' 


(B .21) 



APPENDIX C 


MATERIAL PROPERTIES 

The material properties for graphite-epoxy T300/5208 are given by Nagar- 
kar and Herakovich [25] as 


Elastic Moduli 

Ejj = 19.2 x 106 p S i t £ 22 = 1*56 x 10® psi, E 33 = 1.56 x 10® psi 


Shear Moduli 

G 23 = °.487 x 106 P si » G i 3 = 0*820 x 10® psi, G 12 = 0.820 x 10® psi 


Possion's Ratios 

v 23 = 0.490, v 13 = 0.238, = 0.238 


X T = 219.5 x 10 3 psi, 
X c = -246.0 x 10 3 psi, 
3 23 = g * 8 x ^ psi » 


Strength Parameters 
Yj = 6.35 x 10 3 psi, 

Y c = -23.85 x 10 3 psi , 
S 13 = 12.6 x 10 3 psi. 


Zj = 6.35 x 10 3 psi 
Zq - -23.85 x 10 3 psi 
Sj 2 = 12.6 x 10 3 psi 


The material properties for the isotropic problems considered were 
chosen to be 

E = 30 x 10® psi, v = 0.20, o u = 50 x 10 3 psi 
where is the ultimate strength of the material. 
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FINITE ELEMENT MESHES 
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APPENDIX E 


ON THE RELATIONSHIP BETWEEN FREE-EDGE STRESSES AND 
THE DIRECTION OF CRACK EXTENSION IN ANGLE-PLY LAMINATES 

It was originally thought that the free-edge stress state could be 
used to predetermine the direction of crack extension. 

Herakovich [7, 8], for example, has demonstrated that the crack 

types present on the free-edge of tension specimens seem to exhibit 

distinct crack growth patterns direction of extension for different 

angle-ply laminates. To account for this effect it was thought that the 

orientation of the principal stress plane on the free edge would follow 

along the same path as the crack plane when traced through the laminate 

thickness. To test this premise, the orientation of the crack 

plane, <j> , through the laminate thickness was compared with the through 

the thickness variation of the principal stress plane. A generalized 

plane strain solution of the front face. Fig. (3.6a), was obtained using 

the same geometry, loading and boundary conditions as given in Section 

3.3. Only the stresses in the plane of the free edge were considered, 

i.e., o ,o and x , in Fig. (2.10). It was assumed that the thermal 
xx zz xz 

curing stresses were balanced by the hygroscopic stresses in the lami- 
nate so that only mechanical loading be considered. The material con- 
sidered and experimental crack plane orientation was that of Herakovich 
[7] (the material was the same T300/5208 graphite-epo^y of this report-- 
with properties listed in Appendix C). The orientation of the principal 
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plane, <J> , was calculated through 


* = 270° - |tan _1 [ : 2 xz 


a -a 
xx zz 


(E.l) 


The results for four laminate configurations are presented in Fig. 
(E.la-d). The multi-valued angles, 4>p * at z/H = 0.5 represent the 
effect of a crack turning into a delamination at that point. Fig. 
(E.la-d) indicate that only in the [902/02]$ case did the theory agree 
with the experiment. Hence, it was concluded that the direction of 
crack extension could not be predicted through the crack free stress 
state but that an actual crack must be introduced. 



Figure E.1 Principal Plane Orientation vs. Plane 
Of Crack Extension For Angle-Ply Laminates 








APPENDIX F 


COMPOSITE LAMINATE FRACTURE MECHANICS (CLFM2D) 
INPUT DATA SEQUENCE 


Cards 1 & 2 : title cards (20A4) 

Column Contents 

1-80 ITITLE (I,J) Title 


Card 3 : Control card 

Column 

1-5 NPROB 


6-10 NEM 
11-15 NODS 
16-20 NANG 


(11IS) 

Contents 

Problem type (1 = generalized plane 
strain, 2 = skewed plane stress, 

3 = orthotropic plane stress, 4 = iso 
tropic plane strain) 

Number of elements in mesh 
Number of nodes in mesh 
Number of different angles - must be 
> 1 


21-25 NSDF 
26-30 NSBF 
31-35 NEXX 


Number of specified degrees of freedom 
Number of specified forces (tractions) 
Number of different E xx - must be > 1. 
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36-40 NFX 


41-45 NPLOT 


46-50 NT1 


51-55 NT2 


56-60 NCHECK 


Fracture mechanics key 
NFX = 0, no fracture analysis 
NFX f 0, fracture mechanics problem 
Plot key 

NPLOT = 0, no plot 
NPLOT = 1, undeformed plot only 
NPLOT = 2, deformed plot only 
NPLOT = 3, both deformed and 
undeformed plots 
Dump file for displacements 
NT1 = 0, no dump 
NT1 i 0, dumps displacements to 
unit NT1 

Dump file for strains, stresses and 
strain energies 
NT2 = 0, no dump 
NT2 f 0, dumps to unit NT2 
Data check option (NCHECK f 0 - data 
check only) 


Card 4: 

Scale factors 

Column 


1-10 

SCAY 

11-20 

SCAZ 


(2D10.5) 

Contents 
Y-scale factor 
Z-scale factor 
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Card 5: 
Column 

Printer 

control card (615) 

Contents 

1-5 

KEY (1) 

Key for printing element data 

6-10 

KEY (2) 

Key for printing nodal data 

11-15 

KEY (3) 

Key for printing specified displacement 
data 

16-20 

KEY (4) 

Key for printing specified force data 

21-25 

KEY (5) 

Key for printing nodal displacements 

26-30 

KEY (6) 

Key for printing strains, stresses and 
strain energies 

Note: If KEY 

(I) ^ 0, print 

Card 6: 
Column 

Plotter 

scale factors (4D10.5) - skip if NPLOT = 0 
Contents 

1-10 

PVSCL 

Plot Y-scale factor 

11-20 

PZSCL 

Plot Z-scale factor 

21-30 

VMAX 

Maximum Y-displacement 

31-40 

WMAX 

Maximum Z-di splacement 

Card 7: 
Column 

Material 

property card 1 (6D10.5) 
Contents 

1-10 

PROP (1) 

E 11 

11-20 

PROP (2) 

e 22 

21-30 

PROP (3) 

e 33 

31-40 

PROP (4) 

g 23 


Contents 
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41-50 

PROP (5) 

G 13 

51-60 

PROP (6) 

G 12 

Card 8: 

Material property 

card 2 (6D10.5) 

Column 


Contents 

1-10 

PROP (7) 

v 23 

11-20 

PROP (8) 

v 13 

21-30 

PROP (9) 

°12 

31-40 

PROP (10) 

“ll 

41-50 

PROP (11) 

a 22 

51-60 

PROP (11) 

a 33 

Card 9: 

Angles (8010. 5) 


Column 


Contents 

1-10 

ANG (1) 

Angle No. 1 (in 

11-20 

• 

ANG (2) 

Angle No. 2 (in 

• 

• 

• 

ANG (NANG) 

• 

• 

Angle No. "NANG 


degrees) 

degrees) 


(in degrees) 
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Card 10: Element data {5x, 1415) repeat "NEM" times 

Column 


Contents 

1-5 


Blank (can put in element numbers) 

6-10 

NOD (1,1) 

Node No. 1 of element I* 

11-15 

NOD (1,2) 

Node No. 2 of element I* 

16-20 

NOD (1,3) 

Node No. 3 of element I* 

21-25 

NOD (1,4) 

Node No. 4 of element I* 

26-30 

IANG(I) 

Angle number of element 1^ 

31-35 

IEPS (1,1) 

0 

e xx number at local node 1 

36-40 

IEPS (1,2) 

e xx number at local node 2 

41-45 

IEPS (1,3) 

e xx number at local node 3^ 

46-50 

IEPS (1,4) 

e xx number at local node 4 

51-55 

ISTRS (1,1) 

Stress output location l 4 

56-60 

ISTRS (1,2) 

Stress output location 2^ 

61-65 

ISTRS (1,3) 

Stress output location 3^ 

66-70 

ISTRS (1,4) 

Stress output location 4^ 

71-75 

ISTRS (1,5) 

Stress output location 5 4 

1 . 

Refer to Fig. 

(F.l) for node numbering sequence. 

2. 

Refer to Fig. 

(F.2) for angle orientation. 

3. 

e xx numbers correspond to those of Card 16. 

4. 

Stress output 

locations shown in Fig. (F.l). 


Figure F.1 Local Node Numbers And Stress- 
Strain Output Locations 
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original PAGH \3 


Card 11 

: Nodal 

data (5x,2D10.5) - repeat "NODS" times. 

Column 


Contents 

1-5 


Blank (can put in node numbers) 

6-15 

Y(I) 

Y coordinate of node I 

16-25 

Z(I) 

Z coordinate of node I 


Card 12 : Specified displacement data (115, DIO. 5) - repeat "NSDF" times 


and skip if NSDF = 0. 


Column 

Contents 

1-5 ND 

Node number 


6-10 

I DR 

Direction (1 = x, 2 = y, 3 = z) 

11-20 

UBDF(I) 

Specified displacement value 

Card 13 

: Specified force 

data (215, DIO. 5) - repeat "NSBF" times 


and skip i f NSBF 

= 0 . 

Column 


Contents 

1-5 

ND 

Node number 

6-10 

IDR 

Direction (1 = x, 2 = y, and 3 = z) 

11-20 

UVSF(I) 

Specified force value 


Card 15 : Temperature data (1D10.5) 

Column Contents 


1-10 TEMP 


Temperature change 
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Card 16: Normal (x-di rection) strain data (8D10.5) 


Column 


Contents 

1-10 

EPSX(l) 

e value No. 1 

XX 

11-20 

EPSX ( 2 ) 

e value No. 2 

XX 


• 

• 


• 

• 


• 

• 


EPSX(NEXX) 

e value No. "NEXX" 

XX 

*if NFX 

= 0, end of data 


Card 17 

: Fracture control 

card (615,2010.5) 

Column 


Contents 

1-5 

NODE 

Initial crack tip node number 1 

6-10 

N0DE1 

Secondary node No. I 1 

11-15 

N0DE2 

Secondary node No. 2 1 

16-20 

NSE 

Number of stop elements (must be > 1) 

21-25 

NSN 

Number of skip nodes (must be > 1) 

26-30 

NCRT 

Fracture criteria 
NCRT = 1 Griffith 
NCRT = 2 Sih-strain energy density 
NCRT = 3 3-D point stress 
NCRT = 4 2-D point stress 

31-40 

ANOT 

Initial crack length 

41-50 

TMAX 

Maximum run time (CPU-seconds) 


1. Refer to Fig. (F.3) for crack nodes description. 


129 


ORiGifcAi. i-A6>i & 

OF POOR QUALITY 


NODE 1 NODE2 



Figure F.3 Nodes Defining Crack Tip 
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Card 18 : Stop element data* (1615) 


Column 



Contents 

1-5 

NSTOP 

(1) 

Stop element No. 1 

6-10 

NSTOP 

(2) 

Stop element No. 2 


NSTOP (NSE) Stop element No. ”NSE" 

Card 19 : Skip node data^ (1615) 


Column 



Contents 

1-5 

NSKIP 

(1) 

Skip No. 1 

6-10 

NSKIP 

(2) 

Skip node No. 2 


NSKIP (NSN) Skip node No. "NSN" 

1. Stop elements terminate the solution when they are reached - 
used to prevent tear through. 

2. Skip nodes are used to eliminate the designated nodes as 
crack growth possibilities - used to prevent tear through 
and the crack from growing back on itself. 
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Card 20 

: Strength properties 

(3D10.5) 

Column 


Contents 

1-10 

STRENG (1) 

X T 

11-20 

STRENG (2) 

y t 

21-30 

STRENG (3) 

Zt 


end of data 


